function  stdp_summary()
%
% generate some summary graphs for spike-timing dependent potentioniation
%


global CONTROL

global TC TC2

sf = getmainselection;
if(sf <=  0) 
    return;
end;

pflag = getplotflag;
QueMessage('STDP-Summary', 1); % clear the que

% use the comparison plot figure
hf = findobj('Tag', 'Summary');
if(isempty(hf))
    hf = figure('Tag', 'Summary', 'Name', 'Summary Plots', 'NumberTitle', 'off');
end;
figure(hf);

clf;
QueMessage('STDP', 1); % clear the que

prompt={'dt1','dt2', 'frequency', 'Vmeas (min)'};
def={'-15', '15', '10', '15 20'};
dlgTitle='Selection parameters for STDP';
lineNo=1;
answer=inputdlg(prompt,dlgTitle,lineNo,def);
dt1 = number_arg(answer{1});
dt2 = number_arg(answer{2});
fr = number_arg(answer{3});
mwin = number_arg(answer{4});
%extract info in a general way...
fprintf(1, '\n\n');
tp=strmatch('timing_base', {CONTROL.protocol});
tc=strmatch('tb_cond', {CONTROL.protocol});
ivp=strmatch('ap-iv2', {CONTROL.protocol});
if(isempty(ivp))
    ivp = strmatch('ap_ivl', {CONTROL.protocol});
end;
if(isempty(tc))
    tc = strmatch('timing_cond', {CONTROL.protocol});
end;
%   mwin = [0.75 0.95];
%   mwin = [0.0 0.2];
fprintf(1, 'File\tdt\tSlope1\tAmpl1\thw\t\tP\tdf\n');
TC=[];
TC2=[];
%   dt1 = -50;
%   dt2 = 50;
if(dt2<dt1)
    t = dt2;
    dt2=dt1;
    dt1=t;
end;
mindt=dt1;
maxdt=dt2;
%  fr = 10;
i = 0;
% 9/6/07: new assumption: you select the first timing base protocol
% and we assume that the next conditioning protocol is the one that will be
% used, along with the following timing base protocol. be very careful!!!

for k = 1:length(sf)
    % slope measures
    fn=CONTROL(sf(k)).filename;
    tb=strmatch(fn, {CONTROL.filename});
    jil=sort(intersect(tb, tp)); % protocols with our file.
    jcl=sort(intersect(tb, tc)); % conditioning protocols
%    jiv = sort(intersect(tb, ivp)); % find the iv protocols
    j0 = sf(k); % ji(1); % this is the first baseline protocol
    jc = jcl(find((jcl-sf(k))>0, 1, 'first')); % j1=ji(2); % this is theconditioning
    j1 = jil(find((jil-sf(k))>0, 1, 'first'));
    j2=[];
    if(length(jil) > 2)
%         j2=ji(3); % could be others sometimes...
%     else
%         j2 = [];
     end;
        
% at this point, never refer to any index except j0 and j2 for the data, and j1 for the conditioning        
    dttest = CONTROL(jc).EPSPKdt.shkdt; % spike-epsp time.
    if(dttest >= mindt && dttest <= maxdt && ~isempty(CONTROL(j0).PSP_TC)) % & fr == floor(CONTROL(jc(1)).EPSPKdt.freq))
        % extract slope information
        i = i + 1;
        dt(i) = dttest;
        epspb = CONTROL(j0).PSP_TC.S_epsp_1; % control period
        epsp1 = CONTROL(j1).PSP_TC.S_epsp_1; % first period...
        if(~isempty(j2))
            epsp2 = CONTROL(j2).PSP_TC.S_epsp_1;  % final period (#2 test)
        else
            epsp2 = [];
        end;
        epspt = [CONTROL(j0).PSP_TC.S_epsp_1 epsp2]; % total period (1 and 2)
        [epspc epspcd] = mean_var(epspb); % control period
        psp = number_arg(CONTROL(jc).psp_time);
        latb(i) =  psp(2) - CONTROL(jc).EPSPKdt.shkl;
        t0 = CONTROL(j0).PSP_TC.zt(1); % start of baseline
   %     tc = CONTROL(jc).PSP_TC.zt(1); % maybe should be measured from
   %     start of conditioning? - then need to store trace times in DFILE
        t1 = CONTROL(j1).PSP_TC.zt(1); % start after conditioning
        nepspb = length(epspb);
        nepsp1 = length(epsp1);
        %mwin1= [floor(npr*mwin(1)):floor(npr*mwin(2))]+1;
        tbase = CONTROL(j0).PSP_TC.zt-t0;
        twin1 = CONTROL(j1).PSP_TC.zt-t0; % from start of baseline... 
        tdel = twin1(1); % delay to first post condition sample
        mwin1 = [1:length(tbase)];
        mwin2 = find((twin1-tdel)/60 >= mwin(1) & (twin1-tdel)/60 <= mwin(2));
        nepsp2 = length(epsp2);
        %mwin2 = [floor(npr2*mwin(1)):floor(npr2*mwin(2))]+1;
        % compute amplitudes...
        vepspb = CONTROL(j0).PSP_TC.V_epsp_1; % baseline over full window
        vepsp1 = CONTROL(j1).PSP_TC.V_epsp_1;
        if(~isempty(j2))
            vepsp2 = CONTROL(j2).PSP_TC.V_epsp_1;
        else
            vepsp2 = [];
        end;
       vepspb = vepspb(find(vepspb ~= 0));
       mwin1 = mwin1(find(vepspb ~= 0));
 %      vepsp1 = vepsp1(find(vepsp1 ~= 0));
 %      mwin2 = mwin2 (find(vepsp1 ~= 0));
       vepspt = [CONTROL(j1).PSP_TC.V_epsp_1 vepsp2];
        [vepspc vepspcd] = mean_var(vepspb);
        npv = length(vepspb);
        % extract the spike parameters
        %
        ivp1 = jc; % get spike from the conditioning protocol;
        hw(i) = NaN; % get the half-width of the spikes from the IV (future: use HW of conditioning spikes...)
%         if(~isempty(ivp1) && ~isempty(CONTROL(ivp1).spike) && check_field(CONTROL(ivp1).spike, 'hwa') ...
%                 && ~isempty(CONTROL(ivp1).spike.hwa))
%             hw(i) = CONTROL(ivp1).spike.hwa + CONTROL(ivp1).spike.hwb;
%         else
%             err = stdp_spike_dt(ivp1); % try to do it anyway
%             if(err == 1)
%                 fprintf(1, 'no spike info at %d? (%s, %s)\n', ivp1, CONTROL(ivp1).filename, CONTROL(ivp1).protocol);
%                 CONTROL(ivp1).spike
%             else
%                 hw(i) = CONTROL(ivp1).spike.hwa + CONTROL(ivp1).spike.hwb;
%             end;
%         end;

        mnp = mwin1; %[floor(np*0.3):np]; % baseline duration - last half of baseline
        mnr = mwin2;
        [m, b, r, p] = linreg(mnp/length(mnp),epspb(mnp)/mean_var(epspb(mnp)));
        [m2, b2, r2, p2] = linreg([1:nepspb]/nepspb, epspb/mean_var(epspb));
        epsplat = 0;
        %epsplat = CONTROL(j0).PSP_TC.O_epsp_1;
        epsplat = epsplat + CONTROL(jc(1)).EPSPKdt.shkl;
        %dep(i) = CONTROL(jc(1)).EPSPKdt.depol;
        if(length(jc) > 1 && ~isempty(CONTROL(jc(2)).EPSPKdt) && ~isnan(CONTROL(jc(2)).EPSPKdt.shkdt))
            dt2(i) = CONTROL(jc(2)).EPSPKdt.shkdt;
            %   dep2(i) = CONTROL(jc(2)).EPSPKdt.depol;
        else
            dt2(i)=dt(i); % use the one from the first period - probably there was no second conditioning.
            %   dep(i) = CONTROL(jc(1)).EPSPKdt.depol;
        end;
        [e0(i) e0_s(i)] = mean_var(epspb(mnp));
        [e1(i) e1_s(i)] = mean_var(epsp1(mnr));
        e0_s(i) = sqrt(e0_s(i));
        e1_s(i) = sqrt(e1_s(i));
        if(~isempty(epsp2))
            [e2(i) e2_s(i)] = mean_var(epsp2(mnr2));
            e2_s(i) = sqrt(e2_s(i));
        end;
 
        
        [v1(i) v1_s(i)] = mean_var(vepspb(mnp));
        v1_s(i) = sqrt(v1_s(i));
        [e0(i) e0_s(i)  v1(i)  v1_s(i)];
       epspstd(i) = e0_s(i)/e0(i);
        vepspstd(i) = v1_s(i)/v1(i); % nomralized stdevs.
        m = mean_var(epspb(mnp));
        [mn_e(i) mn_es(i)] = mean_var(epsp1(mnr)/e0(i)); % mean final EPSP slope
        mn_es(i) = sqrt(mn_es(i));
        [vn_e(i) vn_es(i)] = mean_var(vepsp1(mnr)/v1(i)); % voltage
        vn_es(i) = sqrt(vn_es(i));
        if(~isempty(epsp2))
            m = mean_var(epsp1(mnr));
            mv = mean_var(vepsp1(mnr));
            [mn2_e(i) mn2_es(i)] = mean_var(epsp2(mnr2)/m);
            mn2_es(i) = sqrt(mn2_es(i));
            [vn2_e(i) vn2_es(i)] = mean_var(vepsp2(mnr2)/mv);
            vn2_es(i) = sqrt(vn2_es(i));
        end;
        [t, prob, df] = t_test(epsp1(mnr), epspb(mnp));
        vm(i) = CONTROL(j0).Rmp;
        nspk(i) = CONTROL(jc).EPSPKdt.spkn;
        timebase=[(1/6):(1/6):20];
        if(length(epsp1) < length(timebase))
            epsp1(end+1:length(timebase)) = NaN;
        end;
        [TC2(i,:), TC(i,:)] = block_average(timebase, epsp1/mean_var(epspb(mnp)), 1);
        %    fprintf(1, 'file: %s   dt: %7.2f   e1: %7.3f  nspk: %d\n', CONTROL(sf(i)).filename, ...
        %       dt(i), mn_e(i), nspk(i));
        %    if(abs(m) <= 0.2 | abs(m2) <= 0.2)
        %      fprintf('     File %s satifies stability criteria: m = %7.3f, m2 = %7.3f p = %7.3f\n', ...
        %          CONTROL(j0).filename, m, m2, p);
        %   else
        %       fprintf('     File %s fails stability criteria: m = %7.3f, m2 = %7.3f  p = %7.3f\n', ...
        %         CONTROL(j0).filename, m, m2, p);
        %         nspk(i) = -1; % sign that we failed
        %     end;
        fprintf(1, '%15s\t%7.3f\t%8.3f\t%8.3f\t%8.3f\t%9.4f\t%9.4f\t%8.1f\n', ...
            CONTROL(j0).filename, dt(i), mn_e(i), vn_e(i), hw(i), t, prob, df);
        zt = mztime(0, [j0 j1 j2]);
        [P(i), F(i), dfR(i), dfF(i)] = bears_f_test(zt(1:length(epspb)), epspb', zt((length(epspb)+1):length(epspb)+length(epspt)), epspt');
        %      fprintf('      F Test: P = %8.2f   F = %8.2f   dfR = %8.2f   dfF = %8.2f\n', ...
        %         P(i), F(i), dfR(i), dfF(i));

    end;
end;
fprintf(1, '\n');
msloped=mean(epspstd); % find the average of the standard deviation of the baseline region
mvoltd=mean(vepspstd); % for both voltage and slope.

n1 = find(nspk >=0 & nspk < 2.5);
n2 = find(nspk >= 2.5);
n3 = find(nspk < 0); % find stability failures...

ndtp=find(dt>=0); % all positive dt's
ndtn=find(dt<0);
nq{1} = find(mn_e > 1 & dt > 0);
nq{2} = find(mn_e > 1 & dt < 0);
nq{3} = find(mn_e < 1 & dt < 0);
nq{4} = find(mn_e < 1 & dt > 0);

dtp = find(dt>0);
dtn = find(dt<0);
if(length(dtp > 2) & length (dtn > 2))
    [tn, pn, dfn] = t_test(mn_e(find(dt>0)), mn_e(find(dt<0)));
    fprintf(1, '+ vs - dt: P: %8.3f  t: %8.4f  df: %6.0f', pn, tn, dfn);
end;
symlist= ['g', 'b', 'r', 'k'];

subplot(2,4,1);
%plot(dt(nq1), mn_e(nq1), 'go', 'markersize', 2, 'markerfacecolor', 'g');
hold on;
for i = 1:4
  hp = plot(dt(nq{i}), mn_e(nq{i}), 'linestyle', 'none');
    set(hp, 'marker', 'o','markersize', 2, 'markerfacecolor', symlist(i),  ...
        'markeredgecolor', symlist(i), 'linestyle', 'none');
    he = my_errorbar(dt(nq{i}), mn_e(nq{i}), -mn_es(nq{i}), mn_es(nq{i}), symlist(i));
    set(he(end), 'linestyle', 'none');
end;
% plot(dt(nq2), mn_e(nq2), 'bo', 'markersize', 2, 'markerfacecolor', 'b');
% plot(dt(nq3), mn_e(nq3), 'ro', 'markersize', 2, 'markerfacecolor', 'r');
% plot(dt(nq4), mn_e(nq4), 'ko', 'markersize', 2, 'markerfacecolor', 'k');

plot([mindt maxdt], [1-msloped, 1-msloped], 'k--');
plot([mindt maxdt], [1+msloped, 1+msloped], 'k--');
plot([mindt maxdt], [1 1], 'k-');
plot([0 0], [0.0 1.5], 'k-');
u=get(gca, 'YLim');
set(gca, 'YLim', [0 u(2)]);
title('dt vs EPSP slope');

subplot(2,4,2);
%plot(dt(nq1), vn_e(nq1), 'go', 'markersize', 2, 'markerfacecolor', 'g');
hold on
for i = 1:4
    hp =  plot(dt(nq{i}), vn_e(nq{i}), 'linestyle', 'none');
    set(hp, 'markersize', 2, 'markerfacecolor', symlist(i), 'marker', 'o', ...
        'markeredgecolor', symlist(i), 'linestyle', 'none');
    he = my_errorbar(dt(nq{i}), vn_e(nq{i}), -vn_es(nq{i}), vn_es(nq{i}), symlist(i));
    set(he(end),  'linestyle', 'none');
end;
% plot(dt(nq2), vn_e(nq2), 'bo', 'markersize', 2, 'markerfacecolor', 'b');
% plot(dt(nq3), vn_e(nq3), 'ro', 'markersize', 2, 'markerfacecolor', 'r');
% plot(dt(nq4), vn_e(nq4), 'ko', 'markersize', 2, 'markerfacecolor', 'k');
plot([mindt maxdt], [1-mvoltd, 1-mvoltd], 'k--');
plot([mindt maxdt], [1+mvoltd, 1+mvoltd], 'k--');
plot([mindt maxdt], [1 1], 'k-');
plot([0 0], [0.0 1.5], 'k-');
u=get(gca, 'YLim');
set(gca, 'YLim', [0 u(2)]);
title('dt vs EPSP amplitude');

% subplot(2,4,3);
% plot(dt2(nq1), mn2_e(nq1), 'go', 'markersize', 2, 'markerfacecolor', 'g');
% hold on
% plot(dt2(nq2), mn2_e(nq2), 'bo', 'markersize', 2, 'markerfacecolor', 'b');
% plot(dt2(nq3), mn2_e(nq3), 'ro', 'markersize', 2, 'markerfacecolor', 'r');
% plot(dt2(nq4), mn2_e(nq4), 'ko', 'markersize', 2, 'markerfacecolor', 'k');
% plot([-40 40], [1 1], 'k-');
% plot([0 0], [0.0 1.5], 'k-');
% u=get(gca, 'YLim');
% set(gca, 'YLim', [0 u(2)]);
% title('Slope - 2nd Period');
% 
% subplot(2,4,4);
% plot(dt2(nq1), vn2_e(nq1), 'go', 'markersize', 2, 'markerfacecolor', 'g');
% hold on
% plot(dt2(nq2), vn2_e(nq2), 'bo', 'markersize', 2, 'markerfacecolor', 'b');
% plot(dt2(nq3), vn2_e(nq3), 'ro', 'markersize', 2, 'markerfacecolor', 'r');
% plot(dt2(nq4), vn2_e(nq4), 'ko', 'markersize', 2, 'markerfacecolor', 'k');
% plot([-40 40], [1 1], 'k-');
% plot([0 0], [0.0 1.5], 'k-');
% u=get(gca, 'YLim');
% set(gca, 'YLim', [0 u(2)]);
% title('EPSP Amplitude - 2nd Period');


subplot(2,4,5);
plot(e0(nq{1}), e1(nq{1}), 'go', 'markersize', 2, 'markerfacecolor', 'g');
hold on
plot(e0(nq{2}), e1(nq{2}), 'bo', 'markersize', 2, 'markerfacecolor', 'b');
plot(e0(nq{3}), e1(nq{3}), 'ro', 'markersize', 2, 'markerfacecolor', 'r');
plot(e0(nq{4}), e1(nq{4}), 'ko', 'markersize', 2, 'markerfacecolor', 'k');
u=get(gca, 'YLim');
v = get(gca, 'XLim');
mx = max(u(2), v(2));
set(gca, 'Ylim', [0 mx]);
set(gca, 'Xlim', [0 mx]);
plot([0 mx], [0 mx], 'k');
title('baseline slope vs. final slope');


subplot(2,4,6);
% plot(e1(nq{1}), e2(nq{1}), 'go', 'markersize', 2, 'markerfacecolor', 'g');
% hold on
% plot(e1(nq{2}), e2(nq{2}), 'bo', 'markersize', 2, 'markerfacecolor', 'b');
% plot(e1(nq{3}), e2(nq{3}), 'ro', 'markersize', 2,'markerfacecolor', 'r');
% plot(e1(nq{4}), e2(nq{4}), 'ko', 'markersize', 2,'markerfacecolor', 'k');
% %u=get(gca, 'YLim');
% %v = get(gca, 'XLim');
% %mx = max(u(2), v(2));
% set(gca, 'Ylim', [0 mx]);
% set(gca, 'Xlim', [0 mx]);
% plot([0 mx], [0 mx], 'k');
% title('baseline slope vs. final slope');

plot(hw(nq{1}), mn_e(nq{1}), 'go', 'markersize', 2, 'markerfacecolor', 'g');
hold on
plot(hw(nq{2}), mn_e(nq{2}), 'bo', 'markersize', 2, 'markerfacecolor', 'b');
plot(hw(nq{3}), mn_e(nq{3}), 'ro', 'markersize', 2,'markerfacecolor', 'r');
plot(hw(nq{4}), mn_e(nq{4}), 'ko', 'markersize', 2,'markerfacecolor', 'k');
%u=get(gca, 'YLim');
%v = get(gca, 'XLim');
%mx = max(u(2), v(2));
set(gca, 'Ylim', [0 2.5]);
set(gca, 'Xlim', [0 3]);
xlabel('spike hw (ms)');
plot([0 mx], [1 1], 'k');
title('Spike HW vs. final slope');

subplot(2,4,7)
plot(latb(nq{1}), e1(nq{1}), 'go', 'markersize', 2, 'markerfacecolor', 'g');
hold on;
plot(latb(nq{2}), e1(nq{2}), 'bo', 'markersize', 2, 'markerfacecolor', 'b');
plot(latb(nq{3}), e1(nq{3}), 'ro', 'markersize', 2, 'markerfacecolor', 'r');
plot(latb(nq{4}), e1(nq{4}), 'ko', 'markersize', 2, 'markerfacecolor', 'k');
title('latency vs EPSP slope');

subplot(2,4,8)
plot(e1(ndtp), mn_e(ndtp), 'go', 'markersize', 2, 'markerfacecolor', 'g');
hold on;
plot(e1(ndtn), mn_e(ndtn), 'ro', 'markersize', 2, 'markerfacecolor', 'r');
%   plot(latb(nq{3}), v1(nq{3}), 'ro', 'markersize', 2, 'markerfacecolor', 'r');
%   plot(latb(nq{4}), v1(nq{4}), 'ko', 'markersize', 2, 'markerfacecolor', 'k');
title('Potentiation vs EPSP Amplitude');
u=get(gca, 'XLim');
plot(u, [1 1], 'k');


%figure;
n25=find(dt < -12 & dt > -25);
n10=find(dt < 0 & dt >= -12);
p10=find(dt > 0 & dt <= 12);
p25=find(dt > 12 & dt < 25);

[mn25, vn25]=mean_var(mn_e(n25));
sn25 = sqrt(vn25); % std(mn_e(n25));
[mn10, vn10]=mean_var(mn_e(n10));
sn10 = sqrt(vn10); % std(mn_e(n10));
[mp10, vp10]=mean_var(mn_e(p10));
sp10 = sqrt(vp10); % std(mn_e(p10));
[mp25, vp25]=mean_var(mn_e(p25));
sp25 = sqrt(vp25); % std(mn_e(vp25));
% [mn25 sn25]
% [mn10 sn10]
% [mp10 sp10]
% [mp25 sp25]
if(~isempty(n10) & ~isempty(p10))
    [tn, pn, dfn] = t_test(mn_e(p10), mn_e(n10));
    fprintf(1, '+ vs - dt: P: %8.3f  t: %8.4f  df: %6.0f', pn, tn, dfn);
end;
if(~isempty(p25) & ~isempty(p10))
    [tn, pn, dfn] = t_test(mn_e(p25), mn_e(p10));
    fprintf(1, '+ vs - dt: P: %8.3f  t: %8.4f  df: %6.0f', pn, tn, dfn);
end;


function [z] = mztime(k, l)

global CONTROL

z = [];
wrapflag = 0;
for i = l
    j = k + i;
    zt = CONTROL(j).PSP_TC.zt;
    if(wrapflag)
        zt = zt + 60 * 60 * 24;
    end;
    if(i <= 0)
        ts = zt(1)/60; % zero start time, in minutes.
    else
        ts = 0;
    end;
    x=find(diff(zt)<0);
    if(~isempty(x)) % wrap around......
        zt(x:end) = zt(x:end)+60*60*24; % add a day
        %fprintf('wrapflag set\n');
        wrapflag = 1;
    end;
    zt = zt/60;
    if(i > -1)
        t0 = zt(1) - ts; % convert to minutes...
        zt = zt/60+t0;
        z = [z zt];
    else
        z = zt' - ts;
    end;

end;
